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ABSTRACT 


Refinements to the snow-physics scheme of SSiB are described and evaluated. The 
upgrades include a partial re-design of the conceptual architecture to better simulate the 
diurnal temperature of the snow surface. For a deep snowpack, there are two separate 
prognostic temperature snow layers - the top layer responds to diurnal fluctuations in the 
surface forcing, while the deep layer exhibits a slowly varying response. In addition, the use 
of a very deep soil temperature and a treatment of snow aging with its influence on snow 
density is parameterized and evaluated. The upgraded snow scheme produces better timing 
of snow melt in GSWP-style simulations using ISLSCP Initiative I data for 1987-1988 in the 
Russian Wheat Belt region. 

To simulate more realistic runoff in regions with high orographic variability, additional 
improvements are made to SSiB’s soil hydrology. These improvements include an orography- 
based surface runoff scheme as well as interaction with a water table below SSiB’s three soil 
layers. The addition of these parameterizations further help to simulate more realistic runoff 
and accompanying prognostic soil moisture fields in the GSWP-style simulations. 

In intercomparisons of the performance of the new snow-physics SSiB with its earlier 
versions using an 18-year single-site dataset from Valdai, Russia, the version of SSiB de- 
scribed in this paper again produces the earliest onset of snow melt. Soil moisture and deep 
soil temperatures also compare favorably with observations. 
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1. Introduction and Motivation 


The Simplified Simple Biosphere Model (SSiB) (Xue et al., 1991) is a biophysical 
model of land-atmosphere interactions, which was designed to simulate land-surface pro- 
cesses in numerical models realistically. The interactions are calculated from the funda- 
mental governing equations (Sellers et al., 1986) and provide fluxes of radiation absorption, 
reflection, and emission together with momentum, sensible, and latent heat to general cir- 
culation models (GCMs) and regional models. SSiB has been calibrated for a number of 
biomes using observational data taken from several regions of the world. In these cali- 
brations and/or evaluations, atmospheric data serve as external forcing, while the model 
simulates soil/vegetation temperature(s), soil moisture(s), and surface fluxes that are com- 
pared with observations. Thus far, validation datasets include the Russian soil moisture data 
(Robock et al., 1995; Schlosser et al., 1997; Xue et al., 1997), the HAPEX-Mobilhy data from 
France (Xue, Zeng, and Schlosser, 1996), the Cabauw data from Netherlands (Chen et al, 
1997), the Anglo-Brazilian Amazonian Climate Observation Study (ABRACOS) data (Xue 
et al., 1996), the Sahelian Energy Balance Experiment (SEBEX), and HAPEX-Sahel field 
measurement data from Niger (Xue, 1997). Some evaluations were part of the Project for 
Intercomparison of Land-Surface Parameterization Scheme (PILPS, Henderson-Sellers et al., 
1993), while others were part of the Global Soil Wetness Project (GSWP, Dirmeyer et al., 
1999) which used ISLSCP Initiative I data (Meeson et al., 1995; Sellers et al., 1995). SSiB 
was included in an NCEP land-surface model intercomparison (Chen et al., 1996) using the 
First ISLSCP Field Experiment (FIFE) data. These evaluations have not only helped to 
improve SSiB, but also to understand the mechanisms of land-surface processes in different 


1 



parts of the world experiencing different climatic conditions, vegetation, and soils. In this 
way, the ongoing research, development, and evaluation of SSiB has paved the way for more 
complex interaction studies in a model coupled to a GCM or regional model. SSiB is used 
within the Goddard Earth Observing System (GEOS) II GCM (Takacs et al., 1994) as well 
as the Eta model used by NCEP. Other major institutions using SSiB-based models in 
an offline mode and/or coupled to a parent atmospheric model include COLA, JMA, and 
UCLA. 

SSiB, however, had a deficiency in snow melt timing and meltwater infiltration, first 
noted by Robock et al. (1995). Several methods were undertaken to improve the simulation 
of these processes in SiB-based models (Xue et al., 1997; Sellers et al., 1996). However, 
when using the global ISLSCP Initiative I data within GSWP, significant snow-physics and 
meltwater infiltration deficiencies surfaced again (Mocko and Sud, 1998). In fact, the defi- 
ciencies were not unique to SSiB. Snow-physics testing and development for models has also 
been done in several land-surface schemes. Among these are Verseghy (1991), Douville et al. 
(1995), Yang et al. (1997), Loth and Graf (1998a), Liston and Sturm (1998), Desborough 
and Pitman (1998), Jin et al. (1999), Sun et al. (1999), and Smirnova et al. (2000). A 
summary of the current range of snow-physics packages within land-surface schemes is de- 
tailed in Slater et al. (2000). The current design benefited from several ideas and concepts 
discussed in the above papers and have adapted some of them for use in the model. 

Huge deficiencies in snow melt and meltwater infiltration motivated the development 
of an earlier snow- physics scheme described in Sud and Mocko (1999). In this development, 
the snow layer was separated from the top soil layer. This separation allowed non-reflected 
shortwave energy to be absorbed in the snow, or to be transmitted through and absorbed in 
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the soil, which had the effect of keeping the soil warm and “blanketed” under the snowpack. 
Once the snow melt began, the warmer non-frozen soil was able to infiltrate the meltwater 
before the entire snowmelt occurred. Larger meltwater infiltration was accompanied by larger 
GCM-simulated high northern latitudinal evapotranspiration and precipitation as observed 
in and around regions of strong seasonal snow cover (Mocko et al., 1999). 

The primary motivation for additional improvements was to reduce the remaining 
biases in the onset of snow melt. Diagnostic tests revealed that the snow surface was not 
generating the amplitude of diurnal temperature that is observed in the early spring season to 
initiate mid-day melting. Consequently, snowpack slowly warmed until the entire snowpack 
reached the melting temperature at which time it started to melt precipitously. Thus, as 
compared to observations, there was delayed initiation of snowmelt, but once it started, it 
was relatively sudden. This gave the useful clue that deep snowpack needs a diurnal layer, 
and the only way to achieve this is by separating it from the rest of the pack. Accordingly, 
a two layer snow model was designed. 

Besides the above, some other deficiencies in SSiB’s hydrology and fluxes also became 
evident. For instance, SSiB had no provision for generating surface runoff as a function 
of orography. This motivated the development of an orography /runoff function. Similarly, 
modifications were needed to improve the sub-surface runoff (or baseflow). An empirical, but 
reasonable, interaction with the water table was also added. The removal of soil moisture 
by evaporation or transpiration can now occur from any soil layer, although a weighting 
function is in place to favor the previous formulations. These changes were made to improve 
both the simulated soil moisture and runoff as compared to available observations. Finally, 
skin temperatures were often too high in hot dry regions which motivated an examination 
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and partial solution of this problem. 

The goal of this research and development activity was to further improve the verifi- 
able simulated results from SSiB, especially with regards to the new snow-physics package. 
In addition to producing a realistic diurnal cycle of the snowpack, the intended changes were 
designed to improve the simulation of the snow accumulation and snow melt timing, meltwa- 
ter infiltration, snow surface as well as soil temperatures, and runoff. A detailed description 
of all the improvements to the snow-physics in SSiB is found in Section 2. The simulated 
results from all three versions of SSiB (original - Xue et al . , 1991, Sud and Mocko, 1999, 
and that described in this paper) are compared against each other and observations from 
1987-1988 using the ISLSCP Initiative I data in Section 3. The same three model versions 
are also compared against observations using an 18-year catchment dataset from Valdai, 
Russia in Section 4. Conclusions and a discussion are found in Section 5. 

2. Changes from one snow layer version of SSiB 

This section describes the changes made from the one snow layer version of SSiB (Sud 
and Mocko, 1999 - OSL, hereafter) to the two snow layer version of SSiB (TSL, hereafter). 
The first sub-section describes additional improvements made to the snow model, while the 
following two sub-sections describe other changes related to the model’s hydrology and flux 
calculations. Each improvement was separately evaluated. 

a. Changes to the snow model 

Reduced, yet systematic, delays of snow melt suggested the need for additional mod- 
ifications of OSL (Figure 1). The deficiency was due to the inability of OSL to simulate a 
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realistic diurnal temperature. Additionally, a very deep snowpack in the real world would 
not have a constant temperature over its entire depth, from near the insulated soil/snow 
interface to the surface which can radiate and cool freely. Thus a need for introducing a 
diurnal snow layer was evidenced. 

1) Diurnal snow layer 

In the TSL design, the number of prognostic temperature layers is maintained at 
three. Under a deep snowpack, there is no need for a diurnal soil layer because of the 
insulating effect of the snow and very little solar income. However, a diurnal snow layer 
emerges atop the snowpack, which can melt during the day due to strong solar input and 
can cool significantly at night, much more than the rest of the snow beneath this layer. So 
when the snow accumulates enough to activate the snow model (currently, the “trigger” is 
0.005 m, or 5 mm, water equivalent), the snowpack is divided into two layers which start 
with the same temperatures, one for the bulk layer (T s ) and one for the diurnal layer ( T sn ) 
of snow. In time, these two layers evolve differently through their own energy budget and 
heat exchange between them. The diurnal layer has a fixed depth of Z sn = 0.004 m (4 
mm) water equivalent regardless of the total depth of snow, and the bulk layer is comprised 
of the remainder of the depth of the snowpack. For deep snow, the force-restore layer 
(Deardorff, 1978) is brought up to work with snow while the soil layers are joined into one 
(To). Schematics showing the major design features of OSL and TSL are found in Figure 2. 

Scientific justifications to moving the diurnal layer atop the snowpack are as follows. 
This top layer is able to respond to the diurnal changes of the atmosphere if the depth of 
the diurnal layer is chosen to best respond to these changes. Indeed, the top of the snow can 
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cool radiatively at night, without replenishing all the energy by conduction from the bulk 
layer. Thus the diurnal layer works as an insulator for the bulk layer. The opposite holds for 
warm days, particularly with midday sun in the early spring. These warm days can produce 
snow melt in late winter/early spring, which was lacking in OSL. The prognostic equations 
for TSL are discussed in Appendix A. 

The blended soil layer interacts with the bulk snowpack atop virtually the same way 
as in OSL. Heat is exchanged between the soil and snow layers through an assumed very 
thin air gap between them by radiation. This air gap functions to introduce a time scale of 
two or three hours in the exchange of heat between snow and soil. The incoming shortwave 
energy can be reflected, absorbed in the snow (in either of two snow layers) or transmitted 
through the snow and absorbed in the soil. The re-freezing of snow melt on the soil (if below r 
the freezing temperature) which raises the soil heat, as well as the snow melt by ground 
heat flux are left essentially unchanged from OSL. The correction of invoking the saturation 
pressure over ice (rather than over water) with a snowpack is also unaltered. Also as in OSL, 
the scheme holds the snow fluxes constant for melting conditions; this eliminates fictitious 
warming of snow and its influence on saturation vapor pressure in the implicit backward 
solution. For details on these processes, see Sud and Mocko (1999). 

2) Age effect 

Having experimented with a few documented algorithms, the age effect of the snow 
on snow density following Verseghy (1991) is adopted. The snow density of the bulk snow 
layer is assumed to be constant with depth, and increases exponentially with time from the 
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fresh snow value of p s min = 100 kg m 3 to p s max = 300 kg m 3 , following: 

p s (t + At) = (p s (t)-p s MAx) ex P (-°- 24 ^) ■*" Ps MAX’ (!) 

where p s is the density in kg m ~ 3 , r = 86400 sec, and At is the timestep in sec. The snow 
density can become more than p s max as a resu ^ °f re-freezing of rain or meltwater, while 
fresh snowfall lowers the density via a weighted average of the old and new snow. The density 
of the diurnal snow layer is held at 119.6 kg m ~ 3 so as to generate a constant depth diurnal 
layer, which is approximately the p s value from Eq. 1 after a 12-hour half-day period. 

The snow density affects two main parameters of TSL. The snow thermal conductivity 
(/c s ) is a strong function of the snow density. With some sensitivity experiments (not shown), 
it was found that making the snow conductivity a function of snow density has a large effect 
on the evolution of the snowpack. The form of the relation among the available has relatively 
little effect in test simulations, although Loth and Graf (1998b) did show an effect in a 
different set of experiments. Regardless, the following functional form is adopted in TSL 
due to Yen (1981): 


( 2 ) 

where k s is the thermal conductivity of snow in W m~ 1 K~ l , Kj = 2.22 W m -1 K~ l is the 
thermal conductivity of ice, and p w = 1000 kg m~ 3 is the density of water. 
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3) Fractional snow cover 


The fractional snow cover is also a function of the snow density, as well as of the snow 
depth. The water equivalent snow depth multiplied by the density of water and divided by 
the density of snow gives the depth of the snow. If the water equivalent snow depth is below 
10e -5 m, the snow fraction is arbitrarily set to zero. Above 10e -5 m depth, the snow cover 
fraction has a functional form: 
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where F s is the fractional snow cover, Z at is the total water equivalent snow depth in m, a is 
35, and 0 is 2. Figure 3 shows a few representative curves for fractional snow cover as a func- 
tion of water equivalent snow depth based on snow density. This function in TSL is designed 
to increase the fraction slowly at first, then very rapidly before asymptoting to 1.0. The frac- 
tional snow cover then affects the four components of albedo (direct/diffuse, visible/near-IR) 
and the radiation absorption coefficients within SSiB. These radiation parameters are also 
affected by melting snow. If either the diurnal or bulk snow layer temperature is > 272.16 
K, the snow is considered to be in a melting phase or recently melted, which reduces the 
albedos and coefficients to 60% of their values. Similar assumptions were already used in 
SSiB. For details on how these parameters are calculated, see Xue et al. (1991) and Sellers 
et al. (1986). 
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4) Deep soil temperature 


Another significant change to the calculation of soil temperatures is in the bottom 
heat flux due to heat exchange with the infinite ground beneath, which happens regardless 
of the snow cover. A calculation for heat exchange between the deep soil temperature, T d , 
and a constant “very deep” soil temperature, Tdd, is added. A force-restore on an annual 
time scale determines the heat exchange and time rate of change of T d : 

f (6) 

where T dd for each grid box was obtained by doing a multi-year average of T d . Adding Eq. 6 
to the model has the beneficial effect of improving the amplitude of the annual cycle of deep 
temperature; specifically, it cools the deep (and by extension, the diurnal soil temperature 
- T g ) layer during the summer, and warms the soil under the snowpack during the winter. 
The warmer soil temperature can produce or cause earlier melting and meltwater infiltration. 
This modification is physically based and the correction generally improves the simulation. 

b. Changes to the hydrology model 

1) Orography-based runoff 

SSiB did not parameterize a dependence of overland flow for high orography. The 
need for this is evident, although it is difficult to design such a parameterization. In SSiB, 
the only two ways that surface flow could occur were: 1) soil wetness fractions exceeding 
1.0, or 2) due to a treatment of infiltration rates of precipitation (especially convective vs. 
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large-scale) described in Sato et al. (1989). The function of the distribution equation is: 


PI(x) = (P c a c + P p a p ) exp bx +(P C c c + P p Cp ), (7) 

where P c and P p are the convective and large-scale precipitation rates (respectively) in mm 
during the timestep, and PI(x ) can be interpreted as an amount distribution which is used 
to calculate throughfall and infiltration rates. The coefficients (a, b, c ) for convective and 
large-scale precipitation are recast to better match observed runoff values with simulated. 
The old and new values are shown in Table 1. 

The runoff from melting snow cover, which had been excessively high in original SSiB, 
was systematically lessened in OSL. Additionally, the values of runoff were low compared 
to observations from many river basins in the annual time scale. Therefore, a hill-slope 
orography function is constructed which runs off a fraction of the available water (AW - 
primarily, rainfall and meltwater) before it has a chance to enter the soil. The simplest 
functional form in which the runoff fraction is a function of the soil slope did not show 
much promise. However, with the availability of the GTOPO30 digital elevation model 
(DEM) orography heights at a horizontal grid spacing of 30 arc seconds, the runoff fraction 
is reformulated as a function of the standard deviation of the orography height (a h ) in each 
1 deg. by 1 deg. grid box. The a h was calculated from the height at each grid box in the 
GTOPO30 DEM to the mean of the DEM heights in the 1 deg. box. After a global map 
of a h was constructed, the annual runoff deficit (observed runoff minus modeled runoff from 
OSL) was calculated for nearly 30 river basins around the globe with observed runoff data 
for 1987-1988 and a sufficient number of rain gauges in the basin (after Oki et al., 1999). 
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Next, the water available for runoff (AW - precipitation plus meltwater minus snowfall) was 
calculated for the same basins on an annual basis. These values are shown on a scatter plot, 
with the Oh scaled by o h MAX = 1843.47 m, and the runoff deficit scaled by AW. To reduce 
the potential dominance of a few grid boxes in the basin, the Oh was weighted by the AW 
at each grid box, then averaged over the basin, then divided by the basin-averaged AW. A 
fit is drawn to this scattered data and is shown in Figure 4. The final form of the function 
for orography-aided surface runoff is: 


AR = min 


2.07* 


<7h 

a h MAX 


)• 


0.5 


* AW, 


(8) 


where AW can be the available water at any time during the model timestep, and AR is the 
additional runoff generated. The precipitation contribution to AW is calculated after the 
runoff produced via Eq. 7, and the meltwater contribution to AW does not include re-freeze. 

2) Linear reservoir drainage 

Some improvements are also made to the linear reservoir drainage of SSiB. Liston 
et al. (1994) introduced a calculation for discharge from the bottom layer that accounts for 
the spatial variability of the grid box. This equation has been recast as: 


BF =(^t) ( w *~ w v.m) pz *< 


(9) 


where BF is the sub-surface runoff (or baseflow) in m from soil layer 3 during the timestep, 
W 3 is the soil wetness of layer 3, W w yj. is the wilting point wetness (see Mocko and Sud, 
1998), P is the soil porosity, and Z 3 is the depth of soil layer 3 in m. If W 3 is below 
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baseflow does not occur. 


3) Gravitational drainage 

The gravitational drainage in SSiB is also modified accordingly. The equation essen- 
tially is unchanged, other than the drainage had been the minimum of Sellers et al. (1986, 
Eq. 62) and W 3 PZ 3 /At, and now the minimum has been replaced by {W 3 — Wgrav) PZ 3 / At, 
where: 


Wgrav — 




(10) 


is the wetness at which Z 3 is half-saturated, ip s is the soil moisture potential at saturation 
in m, and B is the Clapp and Hornberger (1978) parameter. If W 3 is below IVgrav, this half 
of the minimum is zero, and gravitational drainage does not occur. 

4) Water table interaction 

An important change to SSiB’s hydrology is the addition of a water table interaction 
below soil layer 3. If W 3 happens to drop below W w ji t , a factor is put into place such that 
half of this loss during the timestep is replaced by the water table. Because the assumption 
had been made that the water table height is at the level of the bottom of soil layer 3 , the 
baseflow and gravitational drainage described above will move water back from W 3 to the 
water table during wet periods. An elementary parameterization of the effects of a rising 
and falling water table is thus introduced to keep the very deep soil from going too dry (wet) 
in a prolonged drought (rainy) period. 
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c. Changes to the calculation of SSiB fluxes 


1) Bare soil evaporation 

Several changes are made to the way SSiB calculates surfaces fluxes and how the soil 
responds to these fluxes. The bare soil evaporation and transpiration rates calculated by 
SSiB have remained the same; however, the depth at which the water is removed from the 
soil has been changed. For bare soil evaporation, the water can now be removed from any of 
the three layers - on the assumption that scattered capillaries of various shapes and forms 
(non-water areas) in the soil allow w r ater from the lower layers to pass through to the surface. 
The extraction from each of the three layers is weighted by the fractional weight for each 
layer, divided by the sum of the three weights. The weights ( WTBi ) for each of the three 
layers are: 


WTBi 
WTB 2 
WTB 3 


= exp 


A_g 

Rv Tg 


*w; 


-B 


= < 10 -H 'Tft-JWA)( sr; ^ 1 r) 


* exp 


ips 9 
Rv Tg 


*Wo 


-B 


(11) 

( 12 ) 


(13) 


where g — 9.81 m sec ~ 2 is the acceleration due to gravity at the earth’s surface, R^ = 461.5 
J K~ x kg~ l is the water vapor gas constant, Zj and Wi is the depth over soil layer i and soil 
wetness of layer i respectively, and T g is the top SSiB soil layer temperature. Although some 
soil water extraction can now be realized from layers 2 and 3, the majority of the weighting 
still takes bare soil evaporation loss from soil layer 1. 
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2) Extraction by transpiration 


The extraction of water from the soil that is used in transpiration is similarly modified 
with a weighted fraction taken for each of the layers containing vegetation roots. Addition- 
ally, the weight for each layer ( WTTi ) is zero if the soil moisture potential of that layer (?/>,) 
is less than a minimum of soil moisture potential (^i ow )- These weights are calculated as: 

*i = W B , (14) 

^low = V^ w ut’ (!5) 

WTTi = - ^iow) *-^ROOT> ( 16 ) 

where Z t ROOT the depth of layer i which contains roots. The parameter V’low a ^ so 
used as a minimum when calculating the stomatal resistance parameter shown in Xue et al. 
(1991, Eq. 9). 

3) Buoyancy velocity scale 

The bulk aerodynamic formula for the sensible heat flux in SSiB is modified, after 
Mahrt and Sun (1995, Eqns. 11 and 12). The new formula includes a “buoyancy velocity 
scale” defined by: 


w B - 


A- Zi Ad 


B 



(17) 


where z { is the boundary layer depth in m, and A 0 B is the difference between the virtual 
potential temperature at the surface (9 V ) and the virtual potential temperature of the mixed 
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layer. Their buoyancy velocity scale is added to SSiB’s wind speed at reference height M as: 


M = M+(2 b H wbC^C^), (18) 

where bn = 0.00025, and Ctn and Cun are SSiB’s neutral heat and neutral momentum 
transfer coefficients, respectively. This is done for unstable cases only before the calculation 
of the aerodynamic resistance between the canopy air space and the reference height. The 
minimum reference height wind speed is also changed from 2.0 m sec -1 to 0.1 m sec -1 . 
This is a physically-based upgrade and it has the benefit of increasing sensible heat flux and 
leading to some cooling of anomalously large skin temperatures over hot dry desert regions 
in SSiB. 

Two other changes are made to the SSiB flux calculations. The first is that SSiB biome 
type 11 (desert) in previous versions did not explicitly calculate the stomatal resistance, 
whereas TSL does. The second is that the maximum stomatal resistance for the canopy 
layer for any biome type (and the constant value for all ground story layers) is changed from 
1.0e5 to l.OelO. These changes allow SSiB a wider range to calculate the stomatal resistance 
on its own (especially for desert) reducing the influence of an arbitrary cutoff. 

3. Results with GSWP Global data 

The upgraded version of the snow-physics scheme with SSiB was integrated globally 
for 1987-1988 using the 1 deg. by 1 deg. ISLSCP Initiative I data after a ten-year soil moisture 
spin-up procedure. The detailed procedure followed is the same as for previous versions of 
SSiB, and is detailed in Dirmeyer et al. (1999). Again, the Russia Wheat Belt region (55 - 
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65 N; 20 - 60 E) is a particular area of focus because of the availability of observations of 
soil moisture and runoff in this area. 

Northern Hemisphere observations of snow cover from satellite are also available 
through the ISLSCP Initiative I dataset. The weekly data is retrieved from NOAA satellites 
and is simply a function of snow or no snow cover at the grid box. By taking an area-averaged 
time series of observed snow cover of this region, and comparing it to similar curves from 
the simulated model data, as shown in Figure 5, one can see that TSL does a better job of 
simulating snow melt timing. In addition, TSL also captures some of the mid-winter snow 
melt in Jan 1988 that previous versions did not. This benefit results from the diurnal layer 
of the snowpack which is able to melt realistically during warm winter episodes. 

Over this time period, there are two useful soil moisture validation datasets. One, 
a portion of which includes the Russia Wheat Belt region, is described by Vinnikov and 
Yeserkepova (1991); the other, in the state of Illinois in the central United States, is described 
by Hollinger and Isard (1994). The time resolution of these data are generally a week to ten 
days; however, the Russian data is not available in winter due to soil freezing. This data 
is compared against the simulated soil moisture data in terms of available water in the top 
meter of the soil, which is simply the total water in the top meter, minus the soil’s wilting 
point value. The simulated vs. observed values are shown in Figure 6. In the top part of the 
figure for the Russian region, the effect of the soil moisture spin-up is clearly seen as TSL 
has wetter soil than the original. TSL also better reproduces the spring rise in soil moisture 
due to meltwater infiltration, followed by drying in the summer from evapotranspiration. 
OSL also produces a spring rise, but it is delayed after the peak in the observations. In the 
bottom part of the figure for Illinois, where there is relatively little meltwater generated and 
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the change in the snow-physics has relatively little effect, the simulated soil moisture from 
TSL, OSL, and original SSiB compare well with observations. All versions capture the wet 
1987 and dry 1988 in Illinois. However, TSL is relatively wetter on the annual time scale in 
better agreement with the observations. In addition, TSL has the wettest late summer and 
fall, showing the beneficial effects of the water table interaction. 

The simulated against observed basin-scale runoff is also a worthy tool to validate 
a land-surface scheme. However, to be truly useful, the simulated runoff should be routed 
through a river routing scheme. The simulated runoff by all three versions of the model were 
separately routed using the TRIP river routing due to Oki and Sud (1998) and used in Oki 
et al. (1999). The TRIP-routed runoff for two river basins compared to observations from 
stream flow is shown in Figure 7. These basins were chosen because they had available river 
gauge observations that were close to the soil moisture regions defined in Figure 6. The top 
part of the figure shows the simulated and observed flow in the Volga River, which is in 
the Russia Wheat Belt region. One can see the spuriously high and late runoff in original 
SSiB. TSL has an earlier spring peak in the flow (from meltwater) than OSL. The bottom 
of the figure is in the Mississippi River Basin in the central United States. The annual flow 
is generally similar to OSL, except the June to October flow is higher as a result of the 
modifications to baseflow. 

Another useful land-surface scheme validation variable is skin temperature. This value 
can be retrieved from satellite, where it is measured over the satellite pixel size, or simulated 
by a model, where it is the spatial average of the vegetation and soil skin temperatures 
weighted by the fractional vegetation cover. A useful satellite dataset for this time period is 
TOYS Pathfinder Path A due to Susskind et al. (1997). The skin temperature differences 
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between the model simulated values and the satellite retrievals are shown in Figure 8. The 
upper right panel clearly shows a cold bias in the Northern Hemisphere in original SSiB. 
This is from snow cover being simulated in the model but not observed. OSL in the lower 
left panel removes most of this error, but the closest values to observations are found in TSL 
shown in the lower right panel. Especially in the Russia Wheat Belt region and northwestern 
Canada, TSL is closer than OSL in the simulated skin temperature comparison. 

Additionally, adding the Mahrt and Sun (1995) buoyancy velocity scale to the surface 
wind speed has the effect of partially mitigating a warm bias in surface temperatures between 
OSL and TSL in hot dry regions. A closer look at this result is shown in Figure 9, where in 
regions of the Sahara, Australia, and the western United States, the skin temperatures have 
cooled by about a degree. Other areas of cooling are found in India and the Middle East in 
northern winter (not shown). This calculation uses arbitrary constants and can be improved 
to produce better results. 

4. Results with VALDAI region data 

All three versions of the SSiB were integrated using an 18-year single-site dataset from 
Valdai, Russia. This dataset, which is described in Vinnikov et al. (1996) and Schlosser et al. 
(1997), is a midlatitude grassland and was used for PILPS Phase 2(d) (Schlosser et al., 2000). 
The years of the simulation are 1966-1983 and the area is noted for its deep snowpack in 
winter, followed by a strong spring melt. The data used are from the Usadievskiy catchment 
at Valdai where the long-term measurements were taken. Before the 18-year simulation 
began, each version’s initial soil moisture was separately spun-up using the 1966 forcing 
data until the soil moisture reached quasi-equilibrium. 
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The annual cycle of snow accumulation and snow melt is shown in Figure 10. The 
observations, which are recorded monthly (and after snow events and frequently during 
snow melt) are shown in circles. As in the GSWP simulations, original SSiB simulates much 
delayed snow melt, and TSL has the earliest snow melt. There is some indication that TSL 
may be melting snow too early, although it is difficult to isolate a systematic bias. A better 
picture of this early melt is shown in Figure 11, which is a time series over the entire 18 years. 
For many spring melt periods, TSL is too early. However, TSL does well for some years in 
simulating the maximum snow depth, certainly no worse than OSL. TSL is also shown to 
have periods of mid- winter melt, doing a particularly nice job of reproducing the winter 1971 
snow depth. The early melt in some winters in Valdai may be sensitive to arbitrarily chosen 
values for T^, p s for the diurnal layer, or K t given before or in Appendix A. Additionally, 
the simulated snow albedo by all three versions is around 0.61, where the observed value 
at this site is around 0.75 (Robock et al., 1995). A too low albedo causes additional solar 
energy absorption in the snow and soil and results in incorrectly early melt. The calculation 
of SSiB’s snow albedos will be a focus of future research. 

TSL does a better job of simulating the annual cycle of soil moisture as shown in 
Figure 12. The spring rise and March peak in soil moisture due to meltwater infiltration 
is captured quite well. The peak from OSL is too late, while original SSiB is too dry all 
year with no spring soil moisture peak. TSL (as well as OSL and original SSiB) is too dry 
in late summer, in part due to too high evapotranspiration (shown later). Still, the annual 
average simulated soil moisture is highest with TSL. The 18-year time series of soil moisture 
is seen in Figure 13. Of particular interest is the simulated soil moisture profiles for the years 
1976-1977. TSL keeps the soil moisture wetter as observed for these two years. Evidently, 
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the spring 1976 rise in soil moisture, well reproduced by TSL, wettens the soil and is kept 
at this wetness for this period. 

The annual averaged runoff at the Valdai location is shown in Figure 14. These 
catchment discharges are not routed through TRIP, as the region is sufficiently small. As 
before, the runoff from original SSiB is too high and too late in the season. OSL has a 
very low spring soil moisture peak. In TSL, there is a strong initial peak in the runoff, 
but this drops off suddenly. The late summer and fall runoff is reproduced well, but the 
annual average of TSL is too low for this region. This may be caused by too high simulated 
evapotr anspir at ion . 

The annual cycle of evapotranspiration is shown in Figure 15. TSL well simulates the 
early spring evapotranspiration, as a result of the earlier snow melt and meltwater infiltration. 
However, the summer evapotranspiration in TSL is too high. This change is a combination 
of higher early summer soil moisture (as observed) and the modifications to SSiB described 
earlier. 

Figure 16 shows the annual cycle of simulated SSiB deep temperature (' T d ) against 
the observed soil temperature at 80 cm below the surface. The original model, which had 
the coupled snow/soil layer, is unrealistically cold during the winter. TSL has the smallest 
annual amplitude in better agreement with observations. This suggests the addition of Tdd 
to an annual cycle of force-restore of the deep soil temperature is a useful modification. 

5. Conclusion and discussions 

Overall, the upgraded snow-physics package, along with other improvements to SSiB, 
has been shown to improve SSiB’s simulation of snow depth, meltwater infiltration, runoff, 
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and soil temperatures. Both in the Russia Wheat Belt region and the Valdai simulation, 
the model’s improved physics lead to a better simulation of temperatures, which results in 
earlier melt and well-simulated soil moistures. 

Some problems still remain, however. For one, the model may now be melting too 
early in some regions. The snow melt timing is intimately related to diurnal snow layer 
properties and its thickness. Values taken for p s , tz t , and albedo, which are kept constant 
everywhere, deserve more attention (e.g., Loth and Graf, 1998b). Also, the values taken for 
Tdd have a large effect on the soil temperature under the snowpack and thus the initiation 
of snow melt. The formulation of fractional snow cover and the effect of melting snow on 
albedo also seem to play an influential role, and better parameterizations for these processes 
need to be instituted. In addition, adding the effects of blowing snow (Liston and Sturm, 

1998) , solar absorption by snow on sloped surfaces, and sub-grid snow cover effects (Liston, 

1999) are expected to affect the simulation. 

Additional multi-year validation datasets are needed to test the model. In addition 
to useful single-site datasets such as Valdai, longer time period global datasets would be 
useful to test the model at many different locations and times. The upcoming GSWP 1.5 
and GSWP 2 (using ISLSCP Initiative II data) projects are expected to lead to further 
opportunities to evaluation and improve SSiB’s snow-physics and hydrology. 

SSiB has also been coupled to a very high resolution soil hydrology model for detailed 
soil moisture profiles. In its current form, 100 soil moisture layers of 5 cm thickness for a 
total depth of 5 m are used within SSiB to improve simulation of the inter-layer exchanges 
of water. Very often, the water table is above this 5 m depth, so the water table interaction 
is explicitly resolved. A soil model such as this may lead to improved parameterizations in 
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the three layer SSiB, and improve the simulation of evapotranspiration and baseflow. 

APPENDIX A 

The prognostic equations for the three levels of the two snow layer model of SSiB 
(TSL) are described in this appendix. The subscripts used in this section refer to the layers 
as follows: sn = thin diurnal skin layer of snowpack, s = bulk snowpack layer, and d = 
blended soil and deep soil under snowpack. The prognostic equation for the diurnal, or top, 
snow layer is: 


C H 


dT sr 

1 at 


— Rn (sn) Hsn AE sn 




(T sn - T,) , 


(Al) 


where C sn is the effective heat capacity, ( sn ) is the net radiation, H sn is the sensible heat 
flux, A E sn is the latent heat flux (where A is the latent heat of vaporization), and r = 86400 
sec. The sensible and latent heat fluxes are the same as Xue et al. (1991, Eqns. A5 and A6), 
other than T gs and are replaced by T sn and e*( sn ). 

The prognostic equation for the bulk snow layer is: 

C s ^ = Rn (,) + Kc (T d - T s ) + (T sn - T t ) , (A2) 

Ot T 


where C s and Rn ( s ) are the effective heat capacity and net radiation of the bulk layer, and 
k c is the snow/soil interface conductivity (given later). 

The prognostic equation for the blended soil and deep soil under the snowpack is: 

^ d ~dt ~ ^ ~ K ° ^ d ~ T*) + 365^ ^ dd ~ I'd) > (A 3 ) 
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where Cd and R^ is the effective heat capacity of the soil and net radiation of the soil. 


The net radiations for all three levels with snowpack are: 


Rn(sn) = SVTj^l — a) [1 — exp(— + LW ± [1 — V c {\ — e c )] 



+a s TX(l-e c )-cj s T A n , 

(A4) 

Rn (5) 

= [SlTjXl - a) - Rn (sn) ] (1 - exp (~K t Z s )) , 

(A5) 

Rn (d) 

— STf^l Cx) R n (sn) ■fi'n (s)> 

(A6) 


where SWj.(l — a) is the sum of all four components (direct/diffuse, visible/near-IR) of the 
shortwave incident on the top of the snowpack and not reflected (a is the albedo), LW^ is the 
longwave incident on the top of the snowpack, Z sn = 0.004 m (4 mm) is the water equivalent 
depth of the diurnal snow layer, Z s is the water equivalent depth of the bulk snow layer, 
a s = 8.76 x 10 8 W m~ 2 K~ A is the Stefan-Boltzmann constant, T c is the canopy temperature 
in K, V c is the fractional cover of the canopy vegetation, and e c is the longwave emissivity. 
The transmittance coefficient for shortwave energy through the snow is K t = 25.0, which 
is the same as in Sud and Mocko (1999, Eq. 1). This value allows 10 % of the incoming 
non-reflected solar to transmit through 10 cm of snow/ice. 

The blended conductivity of heat for the snow/soil interface - which conducts heat 
through the layers to the interface and radiates heat through a small assumed air-gap between 
the layers is: 


K c = 1 / 


"o.5(l + \/365 )_Zi 




+ 


+ 


0.5Z s p w 


4a s (^) 3 *sPs 


(A7) 
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where K d is the thermal conductivity of the soil, and k s is the thermal conductivity of snow 
defined in Eq. 2. 

The above equations are solved in a implicit backward method given in Sellers et al. 
(1986) and explained in more detail for use with snow-physics in Sud and Mocko (1999). 
Other than the introduction of the very deep temperature equation (Eq. 6), the prognostic 
equations for temperatures for snow-free land remain the same in SSiB due to Xue et al. 
(1991). 
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Tables: 

Table 1: Current and original values for the coefficients for the amount distribution of 
precipitation in SSiB. 


Version of SSiB 

a c 

b 

Cc 

dp 

Cp 

New 

5.0 

5.0 

6.737946e-3 

0.0001 

0.9999 

Original 

20.0 

20.0 

0.206e-8 

0.0001 

0.9999 
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Figure 1: Line plot showing area-averaged snow cover fraction (0-1) in the Russian Wheat 
Belt region for 1987-1988. The long dashed line is from OSL, while the intermittent dashed 
line is from original SSiB. The thick solid line is from observations taken from satellite. 
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Further additions to OSL SSiB: 

• Snowpack is divided into two layers - a diurnal layer 

atop a bulk layer 

• Only one ground layer acts as an annual force-restore 

layer atop an infinite layer with fixed soil 
temperature 

• Meltwater is further generated by diurnal snow layer 

during warmer part of the day 

References: 

Sud and Mocko, AMS, GEWEX/C.CI P, 1999 


Figure 2: Schematics showing the major design features of two versions of the SSiB snow 
model. On the left is OSL and on the right is TSL. 




Fractional snow cover as function of snow depth 



Figure 3: Fractional snow cover as a function of water equivalent snow depth for several 
representative snow densities: 1) long dash - 100 kg m 3 , 2) short dash - 200 kg m 3 , 3) solid 
- 300 kg m 3 , and 4) intermittent dash - 400 kg m 3 . 
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Figure 4: Scatter plot and least-squares approximation for annual runoff deficit values by 
standard deviation of orography for 29 basins around the globe for 1987-1988 data. 





Fractional Area with Snow Cover in Russian region 



Figure 5: Line plot showing area-averaged snow cover fraction (0-1) in the Russian Wheat 
Belt region for 1987-1988, same as Figure 1, only with the addition of a short dashed line 
which is from TSL. 
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Figure 6: Time series of available water in mm for (top) the Russia Wheat Belt region and 
(bottom) Illinois, United States for TSL, OSL, and original SSiB and observations (along 
with annual averages). 
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Figure 7: Time series of TRIP-routed stream flow in mm day~ l for (top) the Volga River 
Basin in the Russia Wheat Belt region and (bottom) the Mississippi River Basin in the 
central United States for TSL, OSL, and original SSiB and observations from gauges (along 
with annual averages). Also shown are the basin size, the mean orography height, and the 
standard deviation of orography in the basin. 





Figure 8: Four panels showing differences for May 1987 in global monthly averaged skin 
temperatures in K between simulated and observations taken from satellite. The upper left 
panel is from TSL minus original SSiB, the upper right panel is original SSiB minus the 
satellite observations, the lower left panel is OSL minus the satellite observations, and the 
lower right panel is TSL minus the satellite observations. 
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Water Equiv. snow depth for Valdai 



Figure 10: Annual averaged water equivalent snow depth in mm for the 18-year Valdai, 
Russia simulation. TSL is shown with a small dash, OSL is shown with a large dash, and 
original SSiB is shown with a variable dash. The observations are shown as circles plotted 
at the time of year of the measurement. 




Water Equiv. snow depth for Valdai 



Obs TSL Model OSL Model 

Thick/circles Small dash Large dash 


Figure 11: Same as Figure 10, only now showing an 18-year time series of the water equivalent 
snow depth in mm. 




Available water in top meter of the soil for Valdai 



Figure 12: Same as Figure 10, only for the annual averaged soil moisture in mm in the top 
meter of the soil. The observed values are shown with a thick solid line. 




Available water in top meter of the soil for Valdai 



Figure 13: Same as Figure 12, only now showing an 18-year time series of the soil moisture 
in mm in the top meter of the soil. 
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Total Runoff for Valdai 



Figure 14: Same as Figure 12, only for the total runoff in mm day l . The observed values 
are shown with a thin solid line. 






Obs: Thick/solid TSL Model: Small dash OSL Model: Large dash Onq :nl. dosh 

Ave = 1.24605 Ave = 1.50396 Ave = 1.40948 Ave = 1.15053 


Figure 15: Same as Figure 12, only for the total evapotranspiration in mm day 1 




Deep Soil Temperature for Valdai 



Figure 16: Same as Figure 12, only for the deep soil temperature in C. 





